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Abstract:  In  this  study,  we  present  a  novel  modeling  approach  which  combines 
ordinary  differential  equation  (ODE)  modeling  with  logical  rules  to  simulate  an  archetype 
biochemical  network,  the  human  coagulation  cascade.  The  model  consisted  of  five 
differential  equations  augmented  with  several  logical  rules  describing  regulatory  connections 
between  model  components,  and  unmodeled  interactions  in  the  network.  This  formulation 
was  more  than  an  order  of  magnitude  smaller  than  current  coagulation  models,  because 
many  of  the  mechanistic  details  of  coagulation  were  encoded  as  logical  rules.  We  estimated 
an  ensemble  of  likely  model  parameters  ( N  =  20)  from  in  vitro  extrinsic  coagulation  data 
sets,  with  and  without  inhibitors,  by  minimizing  the  residual  between  model  simulations 
and  experimental  measurements  using  particle  swarm  optimization  (PSO).  Each  parameter 
set  in  our  ensemble  corresponded  to  a  unique  particle  in  the  PSO.  We  then  validated  the 
model  ensemble  using  thrombin  data  sets  that  were  not  used  during  training.  The  ensemble 
predicted  thrombin  trajectories  for  conditions  not  used  for  model  training,  including 
thrombin  generation  for  normal  and  hemophilic  coagulation  in  the  presence  of  platelets 
(a  significant  unmodeled  component).  We  then  used  flux  analysis  to  understand  how  the 
network  operated  in  a  variety  of  conditions,  and  global  sensitivity  analysis  to  identify  which 
parameters  controlled  the  performance  of  the  network.  Taken  together,  the  hybrid  approach 
produced  a  surprisingly  predictive  model  given  its  small  size,  suggesting  the  proposed 
framework  could  also  be  used  to  dynamically  model  other  biochemical  networks,  including 
intracellular  metabolic  networks,  gene  expression  programs  or  potentially  even  cell  free 
metabolic  systems. 
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1.  Introduction 

Developing  mathematical  models  of  biochemical  networks  is  a  significant  facet  of  systems  biology. 
Modeling  approaches  differ  in  their  degree  of  detail,  where  the  choice  of  approach  is  often  determined 
by  prior  knowledge,  or  model  requirements  [1],  Ordinary  differential  equation  (ODE)  models  are 
common  tools  for  modeling  biochemical  systems  because  of  their  ability  to  capture  dynamics  and  encode 
mechanism.  However,  ODE  models  typically  come  with  difficult  (or  sometimes  impossible)  parameter 
identification  problems.  For  example,  Gadkar  et  al.  showed  that  even  with  near-perfect  information,  it 
was  often  impossible  to  identify  all  the  parameters  in  typical  signal  transduction  models  [2].  However, 
it  is  not  clear  whether  we  actually  need  precise  estimates  for  all  model  parameters.  Bailey  suggested 
more  than  a  decade  ago,  that  achieving  qualitative  or  even  quantitative  understanding  of  biological 
systems  should  not  require  complete  structural  and  parametric  knowledge  [3].  Since  Bailey’s  complex 
biology  with  no  parameters  hypothesis,  Sethna  showed  that  model  performance  is  typically  sensitive 
to  only  a  few  parameters,  a  characteristic  seemingly  universal  to  multi-parameter  models  referred  to 
as  sloppiness  [4].  Thus,  reasonable  predictions  may  be  possible,  despite  parametric  uncertainty,  if  a 
few  critical  parameters  are  well-defined.  For  example,  Tasseff  et  al.,  showed  in  a  model  of  Retinoic 
acid  (RA)  induced  differentiation  of  HL-60  cells,  that  correct  predictions  were  possible  even  when  75% 
of  the  parameters  were  known  only  to  an  order  of  magnitude  [5].  Perhaps  more  importantly,  ODE 
models  require  significant  mechanistic  information,  thereby  limiting  their  utility  in  poorly  understood 
systems,  or  conversely  explode  in  size  when  considering  multiple  pathways  or  subsystems.  Toward 
this  challenge,  logical  modeling  is  an  emerging  paradigm  that  encodes  causal  relationships  between 
model  components  using  quasi-mechanistic  non-linear  transfer  functions  [6].  Logical  models  are  highly 
flexible,  and  despite  their  simplicity,  they  have  captured  rich  behaviors  in  a  variety  of  systems  important 
to  human  health  [7-9].  However,  modeling  complex  dynamics  with  logical  models  is  challenging. 
Thus,  there  is  an  unmet  need  for  a  third  approach  which  combines  ODEs  and  logical  models,  where 
ODEs  could  encode  mechanistic  information,  while  missing  or  incomplete  mechanistic  knowledge  can 
be  approximated  using  a  logical  approach. 

In  this  study,  we  developed  a  hybrid  approach  which  combined  ODE  modeling  with  logical  rules  to 
model  a  well  studied  biochemical  network,  the  human  coagulation  system.  Coagulation  is  an  archetype 
proteolytic  cascade  involving  both  positive  and  negative  feedback  [10-12].  Coagulation  is  mediated  by 
a  family  proteases  in  the  circulation,  called  factors  and  a  key  group  of  blood  cells,  called  platelets.  The 
central  process  in  coagulation  is  the  conversion  of  prothrombin  (fll),  an  inactive  coagulation  factor,  to 
the  master  protease  thrombin  (Flla).  Thrombin  generation  involves  three  phases,  initiation,  amplification 
and  termination  [13,14].  Initiation  requires  a  trigger  event,  for  example  vessel  injury,  which  leads  to  the 
activation  of  factor  VII  (FVIIa).  Two  converging  pathways,  the  extrinsic  and  intrinsic  cascades,  then 
process  and  amplify  this  initial  coagulation  signal.  The  extrinsic  cascade  is  generally  believed  to  be 
the  main  mechanism  of  thrombinogenesis  in  the  blood  [15-17].  Initially,  thrombin  is  produced  upon 
cleavage  of  prothrombin  by  fluid  phase  activated  factor  X  (FXa),  which  itself  has  been  activated  by  Tissue 
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Factor/activated  factor  VII  (TF/FVIIa)  [10].  Picomolar  amounts  of  thrombin  then  activate  the  cofactors 
factors  V  and  VIII  (fV  and  fVIII)  and  platelets,  leading  to  the  formation  of  the  tenase  and  prothrombinase 
complexes  on  activated  platelets.  These  complexes  amplify  the  early  coagulation  signal  by  further 
activating  FXa,  and  directly  converting  prothrombin  to  thrombin.  There  are  several  control  points  in  the 
cascade  that  inhibit  thrombin  formation,  and  eventually  terminate  thrombin  generation.  Tissue  Factor 
Pathway  Inhibitor  (TFPI)  inhibits  FXa  formation  catalyzed  by  TF/FVIIa,  while  antithrombin  III  (ATIII) 
neutralizes  several  of  the  proteases  generated  during  coagulation,  including  thrombin.  Thrombin  itself 
also  inadvertently  plays  a  role  in  its  own  inhibition;  thrombin,  through  interaction  with  thrombomodulin, 
protein  C  and  endothelial  cell  protein  C  receptor  (EPCR),  converts  protein  C  to  activated  protein  C 
(APC)  which  attenuates  the  coagulation  response  by  proteolytic  cleavage  of  fV/FVa  and  I'VIII/FVIIIa. 
Termination  occurs  after  either  prothrombin  is  consumed,  or  thrombin  formation  is  neutralized  by 
inhibitors  such  as  APC  or  ATIII. 

Previous  coagulation  models  have  typically  been  formulated  as  systems  of  nonlinear  ordinary 
differential  equations,  using  mass  action  or  more  complex  kinetics,  to  describe  the  rates  of  biochemical 
conversions  [18-22].  Mechanistic  ODE  coagulation  models  from  our  laboratory  [23,24]  were  built  upon 
the  earlier  studies  of  lones  and  Mann  [25],  Hockin  et  al.  [26],  and  later  Butcnas  el  al.  [27]  who  developed 
and  then  subsequently  refined  highly  mechanistic  coagulation  models.  Other  laboratories  have  also 
expanded  upon  Hockin  et  al.  for  example  by  exploring  the  intrinsic  pathway,  the  role  of  stochastic 
fluctuations  in  coagulation  [28],  and  the  dynamics  of  thrombin  mediated  clot  formation  [29]  and 
fibrinolysis  [30].  Other  aspects  of  coagulation  have  also  been  modeled,  such  as  platelet  biochemistry  [31], 
multi-scale  models  of  clot  formation  [32,33],  and  transport  inside  clots  [34].  However,  these  previous 
studies  were  largely  based  upon  extensive  mechanistic  knowledge.  This  is  possible  because  blood, 
while  enormously  complex,  can  be  systematically  interrogated.  Other  systems,  such  as  intracellular 
signaling  networks,  are  much  more  difficult  to  experimentally  interrogate.  Towards  this  unmet 
need,  we  formulated  a  hybrid  modeling  approach  which  combines  ODEs  and  logical  rules  to  model 
biochemical  processes  for  which  a  complete  mechanistic  understanding  is  missing.  We  tested  this 
approach  by  modeling  the  human  coagulation  cascade.  Others  have  also  constructed  reduced  order 
human  coagulation  models.  Recently,  Papadopoulos  and  co-workers  constructed  a  phenomenological 
mathematical  model  for  thrombin  generation  [35].  Using  four  ordinary  differential  equations  and  six 
parameters,  they  derived  an  expression  describing  the  temporal  evolution  of  thrombin  generation  in  a 
variety  of  cases.  The  reduced  order  Papadopoulos  model  showed  good  agreement  with  experimental 
data,  and  underscored  that  model  reduction  is  possible  even  for  complex  positive  feedback  systems 
like  coagulation.  However,  the  Papadopoulos  model  was  focused  on  thrombin  generation,  and  had  a 
lesser  emphasis  on  the  influence  of  physiological  inhibitors  such  as  ATIII  or  the  protein  C  pathway.  In 
this  study,  we  focused  on  building  a  reduced  order  coagulation  model  that  included  the  physiological 
inhibitors  using  a  hybrid  strategy.  The  hybrid  model  consisted  of  only  five  differential  equations 
augmented  with  several  logical  rules.  Thus,  the  model  was  more  than  an  order  of  magnitude  smaller 
than  comparable  purely  ODE  models  in  the  literature.  We  estimated  the  model  parameters  from  in  vitro 
extrinsic  coagulation  data  sets,  in  the  presence  of  ATIII,  with  and  without  the  protein  C  pathway.  We  then 
compared  the  model  predictions  with  thrombin  data  sets,  for  both  normal  and  hemophilic  coagulation, 
that  were  not  used  for  model  training.  Once  validated,  we  performed  flux  and  sensitivity  analysis  on 
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the  model  to  estimate  which  parameters  were  critical  to  model  performance  in  several  conditions.  The 
reduced  order  hybrid  approach  produced  a  surprisingly  predictive  coagulation  model,  suggesting  this 
framework  could  potentially  be  used  to  model  other  biochemical  networks  important  to  human  health. 


Figure  1.  Schematic  of  the  connectivity  of  the  reduced  order  coagulation  model.  A  trigger 
compound,  e.g.,  TF/FVIIa  initiates  thrombin  production  (Flla)  from  prothrombin  (fll).  Once 
activated,  thrombin  catalyzes  its  own  activation  (amplification  step),  as  well  as  its  own 
inhibition  via  the  conversion  of  protein  C  to  activated  protein  C  (APC).  APC  and  tissue 
factor  pathway  inhibitor  (TFPI)  inhibit  initiation  and  amplification,  while  antithrombin  III 
(ATIII)  directly  inhibits  thrombin.  All  inhibition  steps  and  trigger-induced  initiation  were 
modeled  using  a  rule-based  approach.  Likewise,  the  dependence  of  amplification  on  other 
coagulation  factors  was  also  modeled  using  a  rule-based  approach.  The  abundance  of  the 
highlighted  species  (in  the  dashed  boxes)  was  governed  by  an  ordinary  differential  equation. 
All  other  species  were  assumed  to  be  constant. 


2.  Results 

2.1.  Formulation  of  Reduced  Order  Coagulation  Models 

We  developed  a  reduced  order  extrinsic  coagulation  model  to  test  our  hybrid  modeling  approach 
(Figure  1).  The  core  of  our  model  was  based  upon  the  earlier  work  of  Ismagilov  and  coworkers  [36-39], 
where  we  added  initiation,  factor  dependence,  and  specific  inhibition  terms  to  the  earlier  simplified 
model.  A  trigger  event  initiates  thrombin  formation  (Flla)  from  prothrombin  (fll)  through  a  lumped 
initiation  step.  This  step  loosely  represents  the  initial  activation  of  thrombin  by  activated  FXa.  Once 
activated,  thrombin  catalyzes  its  own  formation  (amplification  step),  and  inhibition  via  the  conversion  of 
protein  C  to  activated  protein  C  (APC).  Antithrombin  III  (ATIII)  inhibits  amplification,  while  APC  and 
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tissue  factor  pathway  inhibitor  (TFPI)  potentially  inhibit  both  initiation  and  amplification.  All  initiation 
and  inhibition  processes,  as  well  as  the  dependence  of  amplification  upon  other  coagulation  factors, 
was  approximated  using  our  rule-based  approach  (Figure  2).  Individual  regulatory  contributions  to  the 
activity  of  pathway  enzymes  were  integrated  into  control  coefficients  (u’s)  using  an  integration  rule 
(min/max).  These  control  coefficients  then  modified  the  rates  of  model  processes  at  each  time  step. 
Hill-like  transfer  functions  0  <  /  (Z)  <  1  quantified  the  contribution  of  components  upon  a  target 
process.  Components  were  either  individual  inhibitor  or  activator  levels  or  some  function  of  levels,  e.g., 
the  product  of  factor  levels.  In  this  study,  Z  corresponded  to  the  abundance  of  individual  inhibitors 
or  activators,  with  the  exception  of  the  dependence  of  amplification  upon  specific  coagulation  factors 
(modeled  as  the  product  of  factors).  When  a  process  was  potentially  sensitive  to  multiple  inputs,  logical 
integration  rules  were  used  to  select  which  transfer  functions  influenced  the  process  at  any  given  time. 
In  our  proof  of  concept  model,  we  used  a  winner  takes  all  strategy;  the  maximum  or  minimum  transfer 
function  was  selected  at  any  given  time  step.  However,  other  integration  rules  are  certainly  possible. 
Taken  together,  while  the  reduced  order  coagulation  model  encodes  significant  biological  complexity,  it 
is  highly  compact  (consisting  of  only  five  differential  equations).  Thus,  it  will  serve  as  an  excellent  proof 
of  principle  example  to  study  the  reduction  of  a  highly  complex  human  subsystem. 

2.2.  Identification  of  Model  Parameters  Using  Particle  Swarm  Optimization 

A  critical  challenge  for  any  dynamic  model  is  the  estimation  of  kinetic  parameters.  We  estimated 
kinetic  and  control  parameters  simultaneously  from  eight  in  vitro  time-series  coagulation  data  sets 
with  and  without  the  protein  C  pathway.  The  residual  between  model  simulations  and  experimental 
measurements  was  minimized  using  particle  swarm  optimization  (PSO).  A  population  of  particles 
( N  =  20)  was  initialized  with  randomized  kinetic  and  control  parameters  and  allowed  to  search  for 
parameter  vectors  that  minimized  the  residual.  However,  not  all  parameters  were  varied  simultaneously. 
We  partitioned  the  parameter  estimation  problem  into  two  subproblems  based  upon  the  biological 
organization  of  the  training  data;  (i)  estimation  of  parameters  associated  with  thrombin  formation  in  the 
absence  of  the  protein  C  pathway  and  (ii)  estimation  of  parameters  associated  with  the  protein  C  pathway. 
Only  those  parameters  associated  with  each  subproblem  were  varied  during  the  optimization  procedure 
for  that  subproblem,  e.g.,  thrombin  parameters  were  not  varied  during  the  protein  C  subproblem.  The 
PSO  procedure  was  run  for  20  generations  for  each  subproblem,  where  each  generation  was  1200 
iterations.  The  best  particle  from  each  generation  was  used  to  generate  the  particle  population  for  the 
next  generation.  We  rotated  the  subproblems,  starting  with  subproblem  1  in  the  first  generation. 

The  experimental  training  data  for  parameter  estimation  was  reproduced  from  the  experiments  of 
Butenas  and  co-workers  [40].  In  these  experiments  thrombin  generation  was  initiated  by  FVIIa-TF 
using  mean  plasma  concentrations  of  coagulation  proteins  and  inhibitors.  To  prepare  FVIIa-TF,  TF 
(0.5  nmol/L)  was  relipidated  into  400  /imol/L  of  phospholipid  vesicles  (PCPS)  by  incubation  in 
20  mmol/L  HEPES,  150  mmol/L  NaCl,  and  2  mmol/L  CaCl2  pH  7.4  (HBS/Ca2+)  for  30  min  at  37  °C. 
The  relipidated  TF  was  incubated  with  10  pmol/L  factor  Vila  for  20  min  to  allow  the  formation  of 
FVIIa-TF.  Factors  V,  VIII  and  thrombomodulin  (Tm)  (when  protein  C  activation  is  required)  were  added 
to  FVIIa-TF  complex.  Thrombin  generation  was  then  initiated  by  adding  equal  volumes  of  this  mixture 
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with  a  mixture  containing  prothrombin,  factor  IX  and  factor  X,  TFPI,  AT- III  and  protein  C  (added  when 
required),  protein  S  (added  when  required)  and  factor  XI  (added  when  required).  In  the  training  data, 
5  pmol/L  FVIIa-TF  was  used  along  with  200  //mol/L  of  phospholipid  vesicles  (PCPS)  to  initiate 
thrombin  generation.  All  other  the  coagulation  factors  and  inhibitors  i.e.,  factors  X,  IX,  V,  and 
VIII,  prothrombin,  TFPI  and  AT-III,  protein  C  and  protein  S  (when  applicable)  were  at  their  mean 
plasma  concentrations. 
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Figure  2.  Schematic  of  rule -based  effective  control  laws.  Traditional  enzyme  kinetic 
expressions,  e.g.,  Michaelis-Menten  or  multiple  saturation  kinetics  are  multiplied  by  an 
enzyme  activity  control  variable  0  <  Vj  <  1.  Control  variables  are  functions  of 
many  possible  regulatory  factors  encoded  by  arbitrary  transfer  functions  of  the  form 
0  <  fj  (Z)  <  1.  At  each  simulation  time  step,  the  Vj  variables  are  calculated  by  evaluating 
integration  rules  such  as  the  max  or  min  of  the  set  of  transfer  functions  fi, . . . ,  fn  influencing 
the  activity  of  enzyme  E3. 
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Figure  3.  Reduced  order  coagulation  model  training  simulations.  Reduced  order  coagulation 
model  parameters  were  estimated  using  particle  swarm  optimization  (PSO)  without  the 
protein  C  pathway  as  a  function  of  prothrombin.  Solid  lines  denote  the  simulated  mean  value 
of  the  thrombin  profile  for  N  =  20  independent  particles,  points  denote  experimental  data. 

The  shaded  region  denotes  the  99%  confidence  estimate  of  the  mean  simulated  thrombin 
value  (uncertainty  in  the  model  simulation).  (A-C)  depict  training  data  and  results  for  150%, 
100%  and  50%  of  physiological  prothrombin  levels  in  the  absence  of  protein  C  pathway.  The 
experimental  training  data  was  reproduced  from  the  study  of  Butenas  et  al.  [40].  Thrombin 
generation  in  these  experiments  was  initiated  using  5  pmol/L  FVIIa-TF  in  the  presence  of 
200  //mol/L  of  phospholipid  vesicles  (PCPS).  As  depicted  in  (A-C)  the  prothrombin  levels 
were  at  150%,  100%  and  50%  of  their  physiological  concentration  in  the  absence  of  protein 
C  pathway.  All  factors  and  control  proteins  in  these  experiments  were  at  their  physiological 
concentration  unless  otherwise  denoted. 
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Figure  4.  Reduced  order  coagulation  model  training  simulations.  Reduced  order  coagulation 
model  parameters  were  estimated  using  particle  swarm  optimization  (PSO)  with  the 
protein  C  pathway  as  a  function  of  prothrombin.  Only  APC  pathway  parameters  were 
allowed  to  vary  in  these  simulations  keeping  the  parameters  estimated  without  protein  C 
pathways  constant.  Solid  lines  denote  the  simulated  mean  value  of  the  thrombin  profile  for 
N  -  20  independent  particles,  points  denote  experimental  data.  The  shaded  region  denotes 
the  99%  confidence  estimate  of  the  mean  simulated  thrombin  value  (uncertainty  in  the 
model  simulation).  (A-C)  depict  training  data  and  results  for  150%,  100%  and  50%  of 
physiological  prothrombin  levels  in  the  presence  of  the  protein  C  pathway.  The  experimental 
training  data  was  reproduced  from  the  study  of  Butenas  et  al.  [40,41].  Thrombin  generation 
in  these  experiments  was  initiated  using  5  pmol/L  FVIIa-TF  in  the  presence  of  200  //mol/L 
of  phospholipid  vesicles  (PCPS).  As  depicted  in  (A-C)  the  prothrombin  levels  were  at  150%, 
100%  and  50%  of  their  physiological  concentration  in  the  presence  of  protein  C  pathway.  All 
factors  and  control  proteins  in  these  experiments  were  at  their  physiological  concentration 
unless  otherwise  denoted. 


The  reduced  order  coagulation  model  captured  the  role  of  initial  prothrombin  abundance,  and  the 
decay  of  the  thrombin  signal  following  from  ATIII  activity  (Figure  3).  However,  we  systematically 
under-predicted  the  thrombin  peak  and  the  strength  of  ATIII  inhibition  in  this  training  data  set.  On  the 
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other  hand,  with  fixed  thrombin  parameters,  we  captured  peak  thrombin  values  and  the  decay  of  the 
thrombin  signal  (at  least  for  the  150%  fll  case)  in  the  presence  of  both  ATIII  and  the  protein  C  pathway 
(Figure  4).  Lastly,  we  were  unable  to  capture  global  differences  in  initiation  time  across  separate 
data  sets  with  a  single  ensemble  of  model  parameters.  These  differences  likely  resulted  from  normal 
experimental  variability.  For  example,  different  thrombin  generation  experiments  within  the  training 
data  (at  the  same  physiological  factor  levels)  had  significantly  different  initiation  times  (data  not  shown). 
However,  the  inability  to  globally  capture  initiation  time  also  highlighted  a  potential  shortcoming  of 
the  initiation  module  within  the  model.  To  capture  the  variability  in  initiation  time  across  training 
data  sets,  we  included  a  constant  time-delay  parameter  (TD)  for  each  data  group.  The  delay  parameter 
was  constant  within  a  data  set,  but  allowed  to  vary  across  training  data  sets.  Introduction  of  the  delay 
parameter  allowed  the  model  to  simulate  multiple  training  data  sets  using  a  single  ensemble  of  model 
parameters.  Taken  together,  the  model  identification  results  suggested  that  our  hybrid  approach  could 
reproduce  a  panel  of  thrombin  generation  data  sets  in  the  neighborhood  of  physiological  factor  and 
inhibitor  concentrations.  However,  it  was  unclear  whether  the  reduced  order  model  could  predict  new 
data,  without  updating  the  model  parameters. 

2.3.  Validation  of  the  Reduced  Order  Coagulation  Model 

We  tested  the  predictive  power  of  the  reduced  order  coagulation  model  with  validation  data  sets 
not  used  during  model  training.  Two  validation  data  sets  were  used,  thrombin  generation  for  various 
prothrombin  and  ATIII  concentrations  with  the  protein  C  pathway,  and  thrombin  generation  in  normal 
versus  hemophilic  plasma  in  the  presence  of  the  protein  C  pathway.  Lastly,  we  compared  the  qualitative 
output  of  the  model  to  rFVIIa  addition  in  the  presence  of  hemophilia.  The  hemophilia  case  was  an 
especially  difficult  test  as  it  was  taken  from  a  different  study  which  used  a  plasma-based  in  vitro  assay 
involving  platelets  instead  of  phospholipid  vesicles  (PCPS).  All  kinetic  and  control  parameters  were 
fixed  for  the  validation  simulations.  The  only  globally  adjustable  parameter  To,  was  fixed  within 
each  validation  data  set  but  allowed  to  vary  between  data  sets.  The  reduced  order  model  predicted 
the  thrombin  generation  profile  for  ratios  of  prothrombin  and  ATIII  in  the  absence  of  the  protein  C 
pathway  (Figure  5).  Simulations  near  the  physiological  range  (fll, ATIII)  =  (100%,  100%)  or  (125%, 
75%)  tracked  the  measured  thrombin  values  (Figure  5B,C).  On  the  other  hand,  predictions  for  factor 
levels  outside  of  the  physiological  range  (fll, ATIII)  =  (50%,  150%)  or  (150%,  50%),  while  qualitatively 
consistent  with  measured  thrombin  values,  did  show  significant  deviation  from  the  measurements 
(Figure  5A,D).  Likewise,  simulations  of  thrombin  generation  in  normal  versus  hemophilia  (missing  both 
fVIII  and  flX)  were  consistent  with  measured  thrombin  values  (Figure  6).  We  modeled  the  dependence 
of  thrombin  amplification  on  factor  levels  using  a  product  rule  (Z  =  fV  x  fX  x  fVIII  x  flX),  which 
was  then  was  integrated  using  a  min  integration  rule  into  the  control  variable  governing  amplification. 
Thus,  in  the  absence  of  fVIII  or  flX,  the  amplification  control  variable  evaluated  to  zero,  and  the  only 
thrombin  produced  was  from  initiation  (Figure  6B).  However,  the  decay  of  the  thrombin  signal  was 
underpredicted  in  the  normal  case  (Figure  6A),  while  the  activated  thrombin  level  was  overpredicted 
in  hemophilia  simulations,  although  thrombin  generation  was  far  less  than  normal  (Figure  6B).  Taken 
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together,  the  reduced  order  model  performed  well  in  the  physiological  range  of  factors,  even  with 
unmodeled  components  such  as  platelet  activation  in  the  hemophilia  data  set. 
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Figure  5.  Reduced  order  coagulation  model  predictions  versus  experimental  data  for  normal 
coagulation.  The  reduced  order  coagulation  model  parameter  estimates  were  tested  against 
data  not  used  during  model  training.  Simulations  of  different  levels  of  prothrombin  and  ATIII 
were  compared  with  experimental  data  in  the  absence  of  the  protein  C  pathway.  Solid  lines 
denote  the  simulated  mean  value  of  the  thrombin  profile  for  N  =  20  independent  particles, 
points  denote  experimental  data.  The  shaded  region  denotes  the  99%  confidence  estimate  of 
the  mean  simulated  thrombin  value  (uncertainty  in  the  model  simulation).  (A-D)  prediction 
results  for  (FII, ATIII):  (50%,  150%),  (100%,  100%),  (125%,  75%)  and  (150%, 50%)  of 
physiological  prothrombin  and  ATIII  levels  in  the  absence  of  the  protein  C  pathway.  The 
experimental  validation  data  was  reproduced  from  the  study  of  Butenas  etal.  [40].  Thrombin 
generation  in  these  experiments  was  initiated  using  5  pmol/L  FVIIa-TF  in  the  presence  of 
200  //mol/L  of  phospholipid  vesicles  (PCPS).  As  depicted  in  (A-D)  the  prothrombin  and 
ATIII  levels  were  at  (50%,  150%),  (100%,  100%),  (125%,  75%)  and  (150%,  50%)  of  then- 
physiological  concentrations  in  the  absence  of  protein  C  pathway.  All  factors  and  control 
proteins  were  at  their  physiological  concentration  unless  otherswise  denoted. 
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Figure  6.  Reduced  order  coagulation  model  predictions  versus  experimental  data  with  and 
without  coagulation  factors  VIII  (fVIII)  and  IX  (FIX).  The  reduced  order  coagulation  model 
parameter  estimates  were  tested  against  data  not  used  during  model  training.  Simulations 
of  normal  thrombin  formation  with  ATIII  and  the  protein  C  pathway  were  compared  with 
thrombin  formation  in  the  absence  of  fVIII  and  fIX.  Solid  lines  denote  the  simulated  mean 
value  of  the  thrombin  profile  for  N  =  20  independent  particles,  points  denote  experimental 
data.  The  shaded  region  denotes  the  99%  confidence  estimate  of  the  mean  simulated 
thrombin  value  (uncertainty  in  the  model  simulation).  (A,B)  prediction  results  for  normal 
thrombin  generation  and  thrombin  generation  in  hemophilia.  All  factors  and  control  proteins 
were  at  their  physiological  concentration  unless  others  noted.  Coagulation  was  initiated 
with  0.2  nmol/L  FVIIa.  The  experimental  validation  data  was  reproduced  from  the  study  of 
Allen  et  al.  [42]. 
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Figure  7.  Reduced  order  coagulation  model  predictions  of  rFVIIa  administration. 
(A)  Simulations  of  thrombin  formation  in  the  presence  of  ATIII  and  the  protein  C  pathway 
were  conducted  for  a  range  of  trigger  values  (lx-200x  nominal)  in  the  absence  of  fVIII  and 
flX;  (B)  Comparison  of  thrombin  generation  for  normal  versus  hemophilia  for  lOx  nominal 
trigger.  Solid  lines  denote  the  simulated  mean  value  of  the  thrombin  profile  for  N  =  20 
independent  particles.  The  peak  thrombin  time  for  normal  coagulation  (t*)  is  less  than  rFVIIa 
induced  coagulation  in  hemophilia  (£**),  while  the  peak  thrombin  value  was  greater  in  normal 
coagulation.  The  shaded  region  denotes  the  99%  confidence  estimate  of  the  mean  thrombin 
value  (uncertainty  in  the  model  simulation).  All  factors  and  control  proteins  were  at  their 
physiological  concentration  unless  others  noted. 


The  model  ensemble  predicted  a  direct  correlation  between  thrombin  generation  and  rFVIIa 
addition  in  hemophilia  (Figure  7).  In  the  current  model,  we  cannot  distinguish  between  different 
initiation  sources,  e.g.,  TF/FVIIa  versus  rFVIIa,  as  we  have  only  a  single  lumped  initiation 
source  (trigger).  Thus,  we  simulated  the  addition  of  rFVIIa  in  hemophilia  by  removing  fVIII 
and  flX  from  the  model,  and  modulating  the  initial  level  of  trigger.  Simulations  with  a 
baseline  level  of  trigger  were  consistent  with  the  previous  hemophilia  simulations,  where  the 
only  thrombin  produced  was  from  initiation  (Figure  7A,  lx  trigger).  However,  as  we  increased 
the  trigger  strength,  the  thrombin  profile  began  to  approximate  normal  coagulation,  showing  a 
pronounced  peak  albeit  with  a  slower  peak  time  (Figure  7B,  t**  >  £*).  Further  increases  in 
trigger  strength  resulted  in  decreased  thrombin  peak  time  and  increased  maximum  thrombin  values 
(Figure  7A,  50x  trigger).  Thus,  for  large  trigger  values  (200x  trigger),  the  hemophilic  thrombin  profile 
approximated  normal  coagulation,  where  peak  thrombin  was  achieved  shortly  after  administration  and 
95%  of  the  thrombin  was  gone  by  20  min  after  initiation.  We  performed  flux  analysis  to  understand 
how  the  reduced  order  coagulation  model  balanced  initiation,  amplification  and  termination  of  thrombin 
generation  for  normal  and  hemophilic  coagulation.  Analysis  of  the  reaction  flux  through  the  reduced 
order  network  for  thrombin  generation  in  normal,  hemophilia  and  rFVIIa-treated  hemophilia  identified 
three  distinct  operational  modes  (Figure  8).  We  calculated  the  flux  through  four  lumped  reactions, 
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initiation,  amplification,  thrombin-induced  APC  generation  and  total  thrombin  inhibition  (including  both 
APC  and  ATIII  action).  Directly  after  the  addition  of  a  trigger  (e.g.,  TF/FVIIa  or  rFVIIa),  the  lumped 
initiation  flux  was  the  largest  for  all  three  cases.  However,  within  a  few  minutes  enough  thrombin  was 
generated  by  the  initiation  mechanism  to  induce  the  amplification  stage.  During  amplification,  thrombin 
catalyzes  its  own  formation  and  inhibition  by  generating  activated  protein  C  (APC),  a  potent  inhibitor 
of  the  coagulation  cascade.  For  normal  coagulation,  amplification  and  thrombin  inhibition  were  the 
dominate  reactions  by  6  min  after  initiation  (Figure  8,  left).  After  10  min,  the  dominate  reaction  had 
shifted  to  thrombin  inhibition  (both  ATIII  and  APC  action).  In  hemophilia  (missing  both  fVIII  and  flX), 
the  amplification  reaction  did  not  occur,  and  thrombin  was  produced  only  by  initiation  (Figure  8,  center). 
Initiation  was  quickly  inhibited  by  APC,  and  the  thrombin  level  stabilized  (eventually  decaying  at  longer 
times  because  of  ATIII  activity).  Lastly,  when  50 x  trigger  was  used  to  induce  thrombin  formation  in 
hemophilia  (absence  of  IVIII/flX),  initiation  mechanisms  dominated  for  up  to  6  min  following  initiation 
(Figure  8,  right).  Similar  to  hemophilia  alone,  no  amplification  occurred  in  the  50 x  trigger +hemophilia 
case,  and  the  rate  of  thrombin  generation  was  extinguished  by  the  combined  action  of  ATIII  and 
APC.  Taken  together,  the  hybrid  modeling  approach  captured  the  transition  between  the  modes  of 
thrombin  generation,  as  well  as  the  role  that  inhibitors  play  in  attenuating  the  thrombin  generation 
rate.  Thus,  the  transfer  function  approach  encoded  the  inhibitory  logic  of  this  cascade  in  the  absence 
of  specific  mechanism. 
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Figure  8.  Reaction  flux  distribution  as  a  function  of  time  for  thrombin  generation  under 
normal  (left),  hemophilia  (center)  and  rFVIIa  treated  hemophilia  (right).  Reaction  flux  was 
calculated  for  each  particle  at  T  =  0, 4,  6,  8, 10, 12, 14  min  after  the  initiation  of  coagulation. 
Reaction  fluxes  were  calculated  for  each  particle  in  the  parameter  ensemble  ( N  =  20).  Blue 
colors  denote  low  flux  values  while  red  colors  denote  high  flux  values. 
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Figure  9.  Global  sensitivity  analysis  of  the  reduced  order  coagulation  model  with  respect 
to  the  model  parameters.  (A)  Sensitivity  analysis  of  the  thrombin  peak  time  for  different 
prothrombin  levels  (150%,  100%  and  50%  of  the  physiological  value)  as  a  function 
of  activated  protein  C;  (B)  Sensitivity  analysis  of  the  thrombin  exposure  for  different 
prothrombin  levels  (150%,  100%  and  50%  of  the  physiological  value)  as  a  function  of 
activated  protein  C.  Points  denote  the  mean  total  sensitivity  value,  while  the  area  around 
each  point  denotes  the  uncertainty  in  the  sensitivity  value.  The  gray  dashed  line  denotes  the 
45°  degree  diagonal,  if  sensitivity  values  are  equal  for  different  conditions  they  will  lie  on  the 
diagonal.  Sensitivity  values  significantly  above  or  below  the  diagonal  indicate  differentially 
important  model  parameters.  The  radius  of  the  shaded  region  around  each  total  sensitivity 
value  was  the  maximum  uncertainty  in  that  value  estimated  by  the  Sobol  method. 


2.4.  Global  Sensitivity  Analysis  of  the  Reduced  Order  Coagulation  Model 

We  conducted  a  global  sensitivity  analysis  to  estimate  which  parameters  controlled  the  performance  of 
the  reduced  order  model.  We  calculated  the  sensitivity  of  the  time  to  maximum  thrombin  (peak  time)  and 
the  thrombin  exposure  (area  under  the  thrombin  curve)  for  different  levels  of  prothrombin,  and  protein 
C  (Figure  9).  Globally,  41%  of  the  parameters  shifted  in  importance  between  the  (fII,PC)  =  (50%,  0%) 
and  (150%,  100%)  cases  for  the  peak  thrombin  time  (Figure  9A).  The  majority  of  these  shifts  involved 
the  interaction  between  increased  prothrombin  and  the  protein  C  pathway,  while  only  5%  were  directly 
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associated  with  increased  prothrombin  alone.  The  rate  constant  for  thrombin  amplification  was  the 
most  important  parameter  controlling  the  peak  thrombin  time.  While  this  parameter  was  differentially 
important  for  different  prothrombin  levels,  and  in  the  presence  or  absence  of  the  activated  protein  C 
pathway,  it  was  consistently  the  most  sensitive  parameter  in  the  model.  The  saturation  constant  governing 
thrombin  amplification  was  the  second  most  important  parameter,  followed  by  the  initiation  control  gain 
parameter.  Other  important  parameters  influencing  the  thrombin  peak  time  included  the  control  gain  for 
activated  protein  C  formation,  and  the  rate  constant  controlling  ATIII  inhibition  of  thrombin  activity. 
On  the  other  hand,  only  27%  of  the  model  parameters  were  differentially  sensitive  between  the 
(fII,PC)  =  (50%,  0%)  and  (150%,  100%)  cases  for  thrombin  exposure  (Figure  9B).  Of  these  parameters, 
all  of  the  shifts  were  associated  with  the  interplay  between  thrombin  formation  and  the  protein  C 
pathway.  The  rate  constant  controlling  ATIII  inhibition  was  the  most  important  parameter  controlling 
the  thrombin  exposure.  While  this  parameter  was  less  important  in  the  presence  of  protein  C  for 
150%  prothrombin  levels,  it  was  significantly  above  all  other  parameters.  Similar  to  the  peak  time, 
for  150%  prothrombin,  the  control  gain  for  activated  protein  C  formation  was  differentially  important 
along  with  the  rate  constant  controlling  amplification.  However,  the  amplification  parameter  was  much 
less  important  for  thrombin  exposure  vs.  peak  time. 

3.  Discussion 

In  this  study,  we  developed  a  reduced  order  model  of  the  human  coagulation  cascade.  We  modeled 
coagulation  because  it  is  well  studied,  has  a  complex  architecture,  and  has  an  abundance  of  experimental 
data  available  for  model  identification  and  validation.  However,  coagulation  was  just  a  proof  of  concept 
test  of  our  approach.  The  proposed  hybrid  framework  could  also  be  used  to  dynamically  model 
other  biochemical  networks,  including  intracellular  metabolic  networks,  gene  expression  programs 
or  potentially  even  cell  free  metabolic  systems.  The  model  consisted  of  five  differential  equations 
augmented  with  several  logical  rules  describing  regulatory  connections  between  model  components 
and  unmodeled  interactions  in  the  network.  We  estimated  model  parameters  from  in  vitro  extrinsic 
coagulation  data  sets,  in  the  presence  of  ATIII,  with  and  without  the  protein  C  pathway.  To  estimate 
parameters,  the  residual  between  model  simulations  and  experimental  measurements  was  minimized 
using  particle  swarm  optimization  (PSO).  However,  not  all  of  the  model  parameters  were  uniquely 
identifiable,  given  the  training  data.  Instead,  we  estimated  an  ensemble  of  likely  parameter  sets  ( N  =  20) 
from  eight  in  vitro  time-series  coagulation  data  sets  with  and  without  the  protein  C  pathway.  Ensemble 
approaches  have  been  used  previously  for  other  signal  transduction  models  [43-47],  and  for  metabolic 
models  [48]  to  estimate  the  impact  of  poorly  constrained  parameter  values  or  poorly  understood  network 
structure  on  simulation  performance.  Thus,  ensemble  approaches  are  common  in  the  dynamic  modeling 
community.  However,  a  unique  feature  of  the  current  study  is  the  direct  connection  between  our  particle 
swarm  approach,  and  the  parameter  ensemble;  each  particle  in  our  swarm  uniquely  corresponded  to 
a  parameter  set  in  our  ensemble.  Thus,  by  constraining  particles  to  operate  in  different  parameter 
regions,  giving  each  particle  a  different  parameter  combination  to  explore,  or  perhaps  even  suppling 
a  different  model  formulation  to  each  particle  we  can  effectively  traverse  through  complex  parameter 
and  model  spaces.  We  validated  the  ensemble  using  thrombin  data  sets  taken  from  multiple  laboratories 
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for  a  variety  of  experimental  conditions  not  used  during  training.  The  ensemble  predicted  thrombin 
trajectories  for  conditions  not  used  for  model  training,  including  thrombin  generation  for  normal  and 
hemophilic  coagulation  in  the  presence  of  platelets  (a  significant  unmodeled  component).  We  then  used 
flux  analysis  to  understand  how  the  network  operated  in  a  variety  of  conditions,  and  global  sensitivity 
analysis  to  identify  which  parameters  controlled  the  performance  of  the  network.  Flux  analysis  showed 
the  logical  rules  formulation  encoded  the  transitions  between  initiation,  amplification  and  termination 
of  thrombin  generation.  Sensitivity  analysis  suggested  that  the  amplification  rate  constant  was  more 
important  to  the  time  to  peak  thrombin,  while  the  ATIII  inhibition  constant  controlled  thrombin  exposure. 
Taken  together,  the  proposed  hybrid  framework  produced  a  surprisingly  predictive  model,  suggesting  this 
approach  could  be  used  to  effectively  model  other  biochemical  networks  important  to  human  health. 

Malfunctions  in  coagulation  can  have  potentially  fatal  consequences.  Aggressive  clotting  involved 
with  Coronary  Artery  Diseases  (CADs),  collectively  accounts  for  38%  of  all  deaths  in  North 
America  [49].  Coagulation  management  during  surgery  can  also  be  challenging,  particularly  with  the 
increase  in  clinical  use  of  antithrombotic  drugs  [50].  Insufficient  coagulation  due  to  genetic  disorders 
such  as  hemophilia  can  also  result  in  recurrent  bleeding.  The  coagulation  factors  VIII  (fVIII)  and  IX 
(fIX)  are  deficient  in  Hemophilia  A  and  B,  respectively  [51-53].  People  with  mild  hemophilia  have 
5%^f0%  of  the  normal  clotting  factor  levels  while  severe  hemophiliacs  have  <1%  [53].  Hemophilia 
can  be  controlled  with  regular  infusions  of  the  deficient  clotting  factors.  However,  clotting  factor 
replacement  sometimes  leads  to  the  formation  of  fVIII  and  fIX  inhibitors  in  vivo  [54].  Alternatively, 
recombinant  factor  Vila  (rFVIIa)  has  been  used  to  treat  bleeding  disorders  [55,56]  including  hemophilia 
with  and  without  factor  VII  I/I  X  inhibitors  [57].  However,  rFVIIa  requires  frequent  administration 
(every  2-3  h),  and  many  questions  remain  about  its  mechanism  of  action,  its  effective  dosage  [54], 
and  its  overall  utility  for  the  treatment  trauma-associated  hemorrhage  [58].  In  this  study,  we  did  not 
model  rFVIIa-induced  coagulation  directly.  Rather,  we  modeled  a  general  trigger  which  initiated  the 
extrinsic  coagulation  cascade.  Since  we  identified  the  model  using  TF/FVIIa,  inherent  to  our  rFVIIa 
simulations  (and  the  rate  constant  governing  initiation)  was  the  presence  of  TF.  However,  even  with  this 
complication,  the  model  generated  potentially  useful  insight  into  the  rFVIIa  mechanism  of  action,  and 
its  possible  shortcomings  especially  for  the  treatment  of  hemophilia.  The  addition  of  rFVIIa  directly 
activated  thrombin  through  the  initiation  pathway.  However,  no  amplification  of  the  thrombin  signal 
occurred  without  fVIII  or  fIX.  Thus,  the  peak  thrombin  signal  was  lower  than  normal  coagulation,  the 
peak  thrombin  time  was  longer,  and  thrombin  generation  was  eventually  inhibited  by  the  combined  action 
of  ATIII  and  the  protein  C  pathway.  However,  as  the  dose  of  rFVIIa  increased,  the  peak  thrombin  time 
decreased  (eventually  saturating  around  200 x nominal  trigger),  and  the  peak  thrombin  value  increased 
such  that  the  thrombin  profile  resembled  normal  coagulation.  Butenas  et  al.  performed  an  extensive 
in  vitro  study  of  rFVIIa-induced  thrombin  generation  under  normal  and  hemophilic  conditions  [59]. 
They  found  qualitatively  similar  trends,  namely  rFVIIa  restored  normal  coagulation  (even  in  the  absence 
of  TF)  for  large  enough  rFVIIa  doses,  although  rFVIIa-induced  coagulation  in  hemophilia  (even  for 
large  rFVIIa  doses)  lagged  the  normal  profile.  These  results  suggest  that  rFVIIa  administration  alone 
might  not  be  able  to  initiate  normal  coagulation  in  recurrent  bleeding,  unless  the  dosage  is  well  above  a 
critical  threshold.  However,  defining  this  threshold,  which  is  likely  patient  specific,  is  difficult  as  there 
is  tremendous  patient  to  patient  variability  even  with  a  normal  coagulation  phenotype  [60]. 


Processes  2015,  3 


194 


The  hybrid  model  simulations  of  thrombin  formation  were  qualitatively  similar  to  other  published 
coagulation  models.  Many  mathematical  models  of  coagulation  dynamics  have  been  built  upon  the 
Hockin-Mann  model  [26].  For  example,  Brummel-Ziedins  and  co-workers  incorporated  the  Protein  C 
(PC)  pathway  into  the  original  Hockin-Mann  network  to  investigate  thrombin  generation  in  cases  of 
familial  PC  deficiency  [61].  Simulations  using  this  model  showed  that  PC  mutations  caused  elevated 
thrombin  levels  without  changing  the  initiation  time  or  the  slope  of  thrombin  generation.  This  trend 
was  qualitatively  captured  by  our  reduced  order  model.  For  example,  we  showed  decreased  peak 
thrombin  concentration  in  the  presence  of  PC  pathway  and  similar  thrombin  generation  slopes,  although 
the  initiation  time  was  different  as  these  were  different  experimental  trials.  Danforth  et  al.  simulated 
normal  thrombin  generation  using  5  pM  tissue  factor  with  all  other  factors  at  their  mean  physiological 
level  [60].  The  initiation  time  in  this  simulation  was  approximately  4.4  min.  When  predicting  the 
normal  thrombin  generation  curve  using  0.2  nM  of  TF/FVIIa,  we  showed  an  initiation  time  between 
4-5  min  in  a  platelet  based  system,  and  an  initiation  time  of  3-4  min  in  a  PCPS  system  with  5  pM 
TF/FVIIa,  although  these  times  were  largely  dictated  by  the  time  delay  parameter  TD.  Thus,  while 
were  not  able  to  explicitly  predict  initiation  time,  we  did  bracket  the  initiation  time  predicted  by  the 
Danforth  study.  Mitrophanov  et  al.  showed,  in  a  study  exploring  the  mode  of  action  of  rfVIIa  [62], 
that  increasing  amounts  of  rFVIIa  accelerated  the  maximal  slope  of  thrombin  generation,  and  both  the 
peak  thrombin  and  initiation  times.  While  we  failed  to  capture  the  effect  of  rFVIIa  on  initiation  time,  we 
correctly  predicted  that  changes  in  the  rFVIIa  concentration  affected  the  maximal  thrombin  slope  and  the 
propagation  phase.  Lastly,  detailed  mechanistic  coagulation  models,  for  example  the  model  by  Luan  and 
co-workers  [24]  or  the  recent  model  by  Mitrophanov  et  al.  [30],  often  contain  hundreds  of  proteins  and 
interactions.  Thus,  it  is  unlikely  that  the  reduced  order  hybrid  model  will  replicate  all  the  rich  behavior 
of  these  other  models.  However,  for  qualitatively  different  cases  such  as  normal  versus  hemophilia,  the 
hybrid  approach  gave  similar  results.  For  example,  akin  to  the  hybrid  model,  Luan  et  al.  modeled  the 
normal  and  hemophilia  data  from  Allen  [42].  However,  unlike  the  previous  Luan  et  al.  study,  we  used 
this  data  for  validation  rather  than  training.  Hybrid  model  simulations  of  the  Allen  et  al.  data  set  were 
surprisingly  consistent  with  the  Luan  et  al.  model.  For  example,  the  amplification  of  a  normal  thrombin 
signal  were  comparable,  and  in  hemophilia  both  correctly  predicted  decreased  thrombin  amplification. 
Taken  together,  the  hybrid  model  performance  was  similar  to  other  full  scale  mechanistic  models  in  the 
literature,  although  we  consistently  failed  to  predict  initiation  times  across  data  sets. 

The  performance  of  the  proof  of  principle  coagulation  model  was  impressive  given  its  limited  size. 
However,  there  are  several  issues  that  could  be  further  explored.  First,  the  prediction  of  initiation 
time  should  be  investigated.  We  were  able  to  estimate  initiation  time  within  a  data  set,  but  unable  to 
predict  initiation  time  across  independent  data  sets.  This  suggested  that  we  should  update  the  initiation 
module  to  distinguish  between  different  triggers,  e.g.,  TF/FVIIa  versus  rFVIIa  alone,  and  to  include  key 
biological  milestones  such  as  FXa  activation  (a  prerequisite  to  thrombin  formation).  Next,  there  are 
several  additional  biological  modules  that  could  be  added  to  the  core  model  presented  here.  First,  we 
could  include  thrombin-induced  platelet  activation  and  the  role  of  activated  platelets  in  amplification.  We 
captured  thrombin  generation  data  in  the  presence  of  platelets,  however,  the  initial  shape  of  the  activation 
curve  and  the  time-scale  of  activation  was  not  always  consistent  with  the  data.  Platelets  are  activated  by 
thrombin  through  the  cleavage  of  the  extracellular  domain  of  protease-activated  receptors  (PARs)  on  the 
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platelet  surface.  Once  activated,  platelets  play  an  important  role  in  amplification,  and  are  key  mediators 
of  the  positive  feedback  driving  amplification.  Thus,  this  biology  is  a  potentially  important  component 
of  an  expanded  model.  We  should  also  add  the  intrinsic  pathway  to  the  model.  The  intrinsic  pathway  is 
triggered  by  contact  activation  of  the  plasma  protease  factor  XI  (fXI)  by  negatively  charged  surfaces  and 
by  thrombin  and  upstream  factors  such  as  activated  plasma  protease  factor  XII  (FXIIa)  [63,64] .  Activated 
platelets  may  also  release  polyphosphate  which  directly  activates  fXII  [65].  Arguably  a  minor  player  in 
acute  bleeding,  contact  activation  could  also  be  important  in  other  wound  healing  contexts.  Finally,  to 
make  the  model  more  clinically  relevant,  we  should  include  the  biochemical  processes  responsible  for 
clot  formation  and  clot  dissolution  (fibrinolysis).  Clot  formation  is  driven  by  thrombin  activity,  while 
fibrinolysis  is  driven  by  plasmin  activity,  a  key  enzyme  that  cleaves  fibrin  (one  of  the  main  materials 
in  a  clot).  Similar  to  coagulation,  fibrinolysis  is  managed  by  several  activating  and  inhibitory  factors 
which  control  the  balance  between  clot  formation  and  dissolution.  Tissue  plasminogen  activator  (t-PA) 
and  urokinase  activate  plasmin,  along  with  contact  pathway  factors  such  as  fXIa.  On  the  other  hand, 
thrombin  activatable  fibrinolysis  inhibitor  (TAFI)  inhibits  the  degradation  of  fibrin  by  plasmin.  Also, 
similar  to  coagulation,  there  is  considerable  fibrinolysis  and  contact  pathway  data  sets  that  can  be  used 
to  train  the  model.  Lastly,  the  choice  of  max/min  integration  rules  or  the  particular  form  of  the  transfer 
functions  could  be  generalized  to  include  other  rule  types  and  functions.  Theoretically,  an  integration 
rule  is  a  function  whose  domain  is  a  set  of  transfer  function  inputs,  and  whose  range  is  v  G  [0,1]. 
Thus,  integration  rules  other  than  max/min  could  be  used,  such  as  the  mean  or  the  product,  assuming 
the  range  of  the  transfer  functions  is  always  /  G  [0, 1].  Alternative  integration  rules  such  as  the  mean 
might  have  different  properties  which  could  influence  model  identification  or  performance.  For  example, 
a  mean  integration  rule  would  be  differentiable,  which  allows  derivative-based  optimization  approaches 
to  be  used.  The  particular  form  of  the  transfer  function  could  also  be  explored.  We  choose  a  Hill-like 
function  because  of  its  prominence  in  the  systems  and  synthetic  biology  community.  However,  the  only 
mathematical  requirement  for  a  transfer  function  is  that  it  map  a  non-negative  continuous  or  categorical 
variable  into  the  range  /  G  [0, 1].  Thus,  many  types  of  transfer  functions  are  possible. 

4.  Materials  and  Methods 

4.1.  Formulation  and  Solution  of  the  Model  Equations 

We  used  ordinary  differential  equations  (ODEs)  to  model  the  time  evolution  of  proteins  (ay)  in  our 
reduced  order  coagulation  model: 


(1) 


where  1Z  denotes  the  number  of  reactions,  A4  denotes  the  number  of  protein  species  in  the  model. 


The  quantity  ry  (x,  e,  k)  denotes  the  rate  of  reaction  j.  Typically,  reaction  j  is  a  non-linear  function  of 
biochemical  species  abundance,  as  well  as  unknown  kinetic  parameters  k  (/C  x  1).  The  quantity  arj 
denotes  the  stoichiometric  coefficient  for  species  i  in  reaction  j.  If  a,3  >  0,  species  i  is  produced  by 
reaction  j.  Conversely,  if  Oij  <  0,  species  i  is  consumed  by  reaction  j,  while  =  0  indicates  species 
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i  is  not  connected  with  reaction  j.  The  system  material  balances  were  subject  to  the  initial  conditions 
x  (t0)  =  xG,  which  were  specified  by  the  experimental  setup. 

Each  reaction  rate  was  written  as  the  product  of  two  terms,  a  kinetic  term  (fj)  and  a  control  term  (vj) 
that  could  depend  upon  many  regulatory  transfer  functions: 

rj  (x,  e,  k)  =  fjVj  (2) 

We  used  multiple  saturation  kinetics  to  model  the  reaction  term  ff 


where  kjiax  denotes  the  maximum  rate  for  reaction  j,  et  denotes  the  enzyme  abundance  which  catalyzes 
reaction  j,  and  KJS  denotes  the  saturation  constant  for  species  s  in  reaction  j.  The  product  in 
Equation  (3)  was  carried  out  over  the  set  of  reactants  for  reaction  j  (denoted  as  raj).  The  control 
term  Vj  depended  upon  the  combination  of  factors  which  influenced  the  activity  of  enzyme  i.  For  each 
enzyme,  we  used  a  rule-based  approach  to  select  from  competing  control  factors  (Figure  2).  If  an  enzyme 
was  activated  by  m  metabolites,  we  modeled  this  activation  as: 

vj  =  max  (/y  (Z) fmj  (Z))  (4) 

where  0  <  f)j  (Z)  <  1  was  a  regulatory  transfer  function  that  calculated  the  influence  of  metabolite  i 
on  the  activity  of  enzyme  j.  Conversely,  if  enzyme  activity  was  inhibited  by  m  metabolites,  we  modeled 
this  inhibition  as: 

Vj  =  1  -  max  (/y  (Z) , . . . ,  fmj  (Z))  (5) 

Lastly,  if  an  enzyme  had  both  m  activating  and  n  inhibitory  factors,  we  modeled  the  control  term  as: 

Vj  —  min  (uj,  dj)  (6) 
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where: 


Uj 


max  (/y  (Z),.. 

■  >  fmj  (2)) 
j+ 

1  -  max  (/y  (Z) , . . . ,  fnj  (Z)) 


(7) 

(8) 


The  quantities  j+  and  j~  denoted  the  sets  of  activating  and  inhibitory  factors  for  enzyme  j.  If 
a  process  has  no  modifying  factors,  we  set  Vj  =  1.  There  are  many  possible  functional  forms  for 
0  <  fij  (Z)  <  1.  However,  in  this  study,  each  individual  regulatory  transfer  function  took  the  form: 


fi  ( Zj ,  kjj)  = 


*5*7 

1  + 


(9) 


where  Z3  denotes  the  abundance  of  the  j  factor  (e.g.,  metabolite  abundance),  and  ktJ  and  r/  are  control 
parameters.  kl3  was  the  species  gain  parameter,  while  r/  was  a  cooperativity  parameter  (similar  to  a 
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Hill  coefficient).  Applying  the  general  framework  to  the  reduced  coagulation  network  resulted  in  five 
ordinary  differential  equations: 


j,  ( PinttPinit  T  ^arnp^arnp ) 


dx  i 

dt 
dx  2 

dt 

dt 
dx  4 

dt 

dxs 

dt 
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(10) 

(11) 

(12) 

(13) 

(14) 


where  x  =  (///,  FI  la ,  PC,  APC,  ATIII)1 .  The  terms  r*u*  in  the  balance  equations  denote  corrected 
kinetic  expressions  for  initiation,  amplification  and  inhibition  processes.  The  rate  of  initiation  rinit  was 
modeled  as: 


rinit  =  Knit  ( trigger ) 


X\ 


K-init,f  1 1  T 


(15) 


where  Knit ,  Kinitjn  are  the  rate  and  saturation  constants  governing  initiation,  respectively.  The  rate  of 
initiation  was  modified  by  vimt,  the  control  parameter  governing  initiation.  Initiation  was  sensitive  to  the 
level  of  trigger  (activator)  and  TFPI  (inhibitor): 


vinit  =  min  (frnit  {TFPI) ,  f+it  {trigger)) 


(16) 


where  the  transfer  functions  /  took  the  form  of  Equation  (9).  The  rate  of  thrombin  amplification  was 
given  by: 


ramp  karnp  {Xl2  ) 


Xi 


K amp, f II  F 


(17) 


where  kamp,  KamPjii  denote  the  rate  and  saturation  constants  governing  amplification,  respectively.  The 
amplification  control  term,  which  modified  amplification  rate,  was  modeled  as  a  combination  of  multiple 
inhibition  terms  and  one  activation  term: 


Vamp  —  min  ( famp  {TFPI)  ,  famp  (#4)  ,  famp  {^ amp ))  (18) 

where  Zamp  =  fV  x  fX  x  fVIII  x  fIX.  Although  f+  {Zamp)  is  an  activating  term,  we  included 
it  in  the  min  integration  rule;  the  factors  in  Zamp  were  essential  for  amplification  (if  any  of  these  factors 
was  missing  the  amplification  reaction  would  not  occur).  Thus,  the  factors  in  Zamp  were  required 
components,  a  classification  that  we  implemented  by  the  min  selection  rule.  The  rate  activated  protein  C 
formation  was  given  by: 

rape  k  APC,  formation  {TM)  — - - - - -  (19) 

formation, PC  T  Aj 

where  Aupc, formation  and  K formation, pc  denote  the  rate  and  saturation  constants  governing  activated 
protein  C  formation,  respectively  and  TM  denotes  the  thrombomodulin  abundance.  We  modeled  the 
control  term  which  governed  APC  formation  as  a  single  thrombin-dependent  activation  term: 


Vapc  =  max  (/+c  (x2)) 


(20) 
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Lastly,  we  included  direct  irreversible  inhibition  of  Flla  by  ATIII: 


(21) 


where  7  was  estimated  to  be  7  =  1.26.  For  ATIII  inhibition  of  Flla,  the  control  variables  Vinh.  atiii 
was  taken  to  be  unity.  The  model  equations  were  encoded  using  the  Python  programming  language  and 
solved  using  the  ODEINT  routine  of  the  SciPy  module  [66].  The  model  files  can  be  downloaded  from 
http  ://www.  varnerlab.  org . 

4.2.  Estimation  of  Model  Parameters  From  Experimental  Data 

Model  parameters  were  estimated  by  minimizing  the  difference  between  simulations  and 
experimental  thrombin  measurements  (squared  residual): 


(22) 


where  Xj  (r)  denotes  the  measured  value  of  species  j  at  time  r,  x3  (r,  k)  denotes  the  simulated  value 
for  species  j  at  time  r,  and  Uj  (r)  denotes  the  experimental  measurement  variance  for  species  j  at  time 
r.  The  outer  summation  is  with  respect  to  time,  while  the  inner  summation  is  with  respect  to  state. 
We  minimized  the  model  residual  using  Particle  swarm  optimization  (PSO)  [67].  PSO  uses  a  swarming 
metaheuristic  to  explore  parameter  spaces.  A  strength  of  PSO  is  its  ability  to  find  the  global  minimum, 
even  in  the  presence  of  potentially  many  local  minima,  by  communicating  the  local  error  landscape 
experienced  by  each  particle  collectively  to  the  swarm.  Thus,  PSO  acts  both  as  a  local  and  a  global 
search  algorithm.  For  each  iteration,  particles  in  the  swarm  compute  their  local  error  by  evaluating  the 
model  equations  using  their  specific  parameter  vector  realization.  From  each  of  these  local  points,  a 
globally  best  error  is  identified.  Both  the  local  and  global  error  are  then  used  to  update  the  parameter 
estimates  of  each  particle  using  the  rules: 


A;  —  9\A.i  +  62^1  ( Ci  —  kj)  +  93r2  ( Q  —  kj) 
k,  =  k*  +  A; 


(23) 

(24) 


where  (#1,  02.  93)  are  adjustable  parameters,  C,  denotes  the  local  best  solution  found  by  particle  i,  and 


Q  denotes  the  best  solution  found  over  the  entire  population  of  particles.  The  quantities  77  and  r2 
denote  uniform  random  vectors  with  the  same  dimension  as  the  number  of  unknown  model  parameters 


(/C  x  1).  In  thus  study,  we  used  (93,  92,  93 )  =  (1.0, 0.05564,  0.02886).  The  quality  of  parameter  estimates 
was  measured  using  goodness  of  fit  (model  residual).  The  particle  swarm  optimization  routine  was 
implemented  in  the  Python  programming  language.  All  plots  were  made  using  the  Matplotlib  module  of 
Python  [68]. 

4.3.  Global  Sensitivity  Analysis  of  Model  Performance 


We  conducted  a  global  sensitivity  analysis,  using  the  variance-based  method  of  Sobol,  to  estimate 
which  parameters  controlled  the  performance  of  the  reduced  order  model  [69].  We  computed  the  total 
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sensitivity  index  of  each  parameter  relative  to  two  performance  objectives,  the  peak  thrombin  time  and 
the  area  under  the  thrombin  curve  (thrombin  exposure).  We  established  the  sampling  bounds  for  each 
parameter  from  the  minimum  and  maximum  value  of  that  parameter  in  the  parameter  set  ensemble.  We 
used  the  sampling  method  of  Saltelli  el  al.  [70]  to  compute  a  family  of  N  (2 d  +  2)  parameter  sets  which 
obeyed  our  parameter  ranges,  where  N  was  the  number  of  trials,  and  d  was  the  number  of  parameters 
in  the  model.  In  our  case,  N  =  10,000  and  d  =  22,  so  the  total  sensitivity  indices  were  computed  from 
460,000  model  evaluations.  The  variance-based  sensitivity  analysis  was  conducted  using  the  SALib 
module  encoded  in  the  Python  programming  language  [71]. 
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